Source code and data found at https://github.com/maRce10/budgie_inbre_stress_vocal_output
Load packages
## add 'developer' to packages to be installed from github
x <- c("remotes", "lubridate", "readxl", "pbapply", "viridis", "ggplot2",
"kableExtra", "knitr", "formatR", "MASS", "sp", "GGally", "brms",
"lme4", "dplyr", "purrr", "forcats", "tidyr", "modelr", "tidybayes",
"cowplot", "ggrepel", "posterior", "ggridges", "maRce10/PhenotypeSpace")
source("~/Dropbox/R_package_testing/sketchy/R/load_packages.R")
load_packages(x)
source("~/Dropbox/R_package_testing/brmsish/R/html_summary.R")
source("~/Dropbox/R_package_testing/brmsish/R/helpers.R")
source("~/Dropbox/R_package_testing/brmsish/R/read_summary.R")
source("~/Dropbox/R_package_testing/brmsish/R/check_rds_models.R")
Functions and global parameters
opts_knit$set(root.dir = "..")
# set evaluation false
opts_chunk$set(fig.width = 10, fig.height = 6, warning = FALSE, message = FALSE,
tidy = TRUE)
read_excel_df <- function(...) data.frame(read_excel(...))
# for reading months in english format
sl <- Sys.setlocale(locale = "en_US.UTF-8")
standard_error <- function(x) sd(x)/sqrt(length(x))
cols <- viridis(10, alpha = 0.7)
col_pointrange <- cols[7]
# read ext sel tab calls
sels <- read.csv("./data/processed/tailored_budgie_calls_sel_tab.csv")
# keep only spectrographic parameters
sels <- sels[, c("sound.files", "selec", "duration", "meanfreq", "sd",
"freq.median", "freq.IQR", "time.IQR", "skew", "kurt", "sp.ent",
"time.ent", "entropy", "meandom", "mindom", "maxdom", "dfrange",
"modindx", "meanpeakf")]
sels$ID <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][1])
sels$month <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][2])
sels$day <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][3])
sels$year <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][4])
sels$date <- paste(sels$day, substr(sels$month, 0, 3), sels$year,
sep = "-")
sels$date <- as.Date(sels$date, format = "%d-%b-%Y")
# acoustic measurements
areas_by_week <- readRDS("./data/processed/acoustic_space_area_by_individual_and_week.RDS")
indiv_ovlp <- readRDS("./data/processed/acoustic_space_density_overlap_to_first_week_by_individual.RDS")
indiv_ovlp$treatment <- factor(indiv_ovlp$treatment, levels = c("Control",
"Medium Stress", "High Stress"))
group_ovlp <- readRDS("./data/processed/acoustic_space_density_overlap_to_group_by_week.RDS")
# group_ovlp <- do.call(rbind, lapply(unique(group_ovlp$ID),
# function(x) { X <- group_ovlp[group_ovlp$ID == x, ]
# X$overlap.to.group <- X$overlap.to.group -
# X$overlap.to.group[X$week ==
# min(as.numeric(as.character(X$week)))] X$distance.to.group <-
# X$distance.to.group - X$distance.to.group[X$week ==
# min(as.numeric(as.character(X$week)))] return(X) }))
group_ovlp$treatment <- factor(group_ovlp$treatment, levels = c("Control",
"Medium Stress", "High Stress"))
Counts per individual
df <- as.data.frame(table(sels$ID))
names(df) <- c("ID", "Sample_size")
df <- df[order(df$Sample_size, decreasing = FALSE), ]
kb <- kable(df, row.names = FALSE)
kb <- kable_styling(kb, bootstrap_options = c("striped", "hover",
"condensed", "responsive"))
print(kb)
|
ID
|
Sample_size
|
|
132YGMM
|
6
|
|
125YGMM
|
12
|
|
160YGHM
|
16
|
|
310YGHM
|
21
|
|
393YGHM
|
25
|
|
303YGHM
|
37
|
|
398WBHM
|
41
|
|
365YGHM
|
44
|
|
399YGLM
|
46
|
|
300YGHM
|
47
|
|
270YGHM
|
64
|
|
407YGHM
|
97
|
|
386WBHM
|
100
|
|
377WWLM
|
108
|
|
367WWNM
|
119
|
|
371YYLM
|
149
|
|
397YGHM
|
175
|
|
378YYLM
|
195
|
|
404WBHM
|
196
|
|
258YGHM
|
223
|
|
389WWLM
|
230
|
|
262YGHM
|
306
|
|
400YGHM
|
360
|
|
388YGLM
|
377
|
|
327YYLM
|
404
|
|
396YBHM
|
455
|
|
403WBHM
|
566
|
|
356WBHM
|
761
|
|
361YGHM
|
767
|
|
402YGHM
|
849
|
|
395WBHM
|
854
|
|
360YGHM
|
900
|
|
390WBHM
|
935
|
|
405YBHM
|
975
|
|
385YBHM
|
1256
|
|
394YBHM
|
1693
|
Physiological parameters
Barplot and effect sizes graph
physio_models <- readRDS("./data/processed/physiological_response_models.RDS")
breath.count <- stack(metadat[, c("D3.Breath.Count", "D7.Breath.Count",
"D14.Breath.Count", "D21.Breath.Count", "D28.Breath.Count")])
weight <- stack(metadat[, c("D3.Bird.Weight..g.", "D7.Bird.Weight..g.",
"D14.Bird.Weight..g.", "D21.Bird.Weight..g.", "D28.Bird.Weight..g.")])
cort.stress <- stack(metadat[, c("D3.CORT.Stress", "D7.CORT.Stress",
"D14.CORT.Stress", "D21.CORT.Stress", "D28.CORT.Stress")])
cort.baseline <- stack(metadat[, c("D3.CORT.Baseline", "D7.CORT.Baseline",
"D14.CORT.Baseline", "D21.CORT.Baseline", "D28.CORT.Baseline")])
stress <- data.frame(metadat[, c("Bird.ID", "Treatment", "Round",
"Treatment.Room")], week = breath.count$ind, breath.count = breath.count$values,
weight = weight$values, cort.stress = cort.stress$values, cort.baseline = cort.baseline$values,
stress.response = cort.stress$values - cort.baseline$values)
# head(stress)
stress$week <- factor(sapply(strsplit(as.character(stress$week), "\\."),
"[[", 1), levels = c("D3", "D7", "D14", "D21", "D28"))
stress$day <- as.numeric(gsub("D", "", as.character(stress$week)))
stress$round <- paste("Round", stress$Round)
stress$treatment <- factor(stress$Treatment, levels = c("Control",
"Medium Stress", "High Stress"))
# remove 5th week
stress <- stress[stress$week != "D28", ]
stress_l <- lapply(stress$Bird.ID, function(x) {
X <- stress[stress$Bird.ID == x, ]
X$breath.count <- X$breath.count - X$breath.count[X$week == "D3"]
X$weight <- X$weight - X$weight[X$week == "D3"]
X$cort.stress <- X$cort.stress - X$cort.stress[X$week == "D3"]
X$cort.baseline <- X$cort.baseline - X$cort.baseline[X$week ==
"D3"]
X$stress.response <- X$stress.response - X$stress.response[X$week ==
"D3"]
return(X)
})
stress <- do.call(rbind, stress_l)
agg_stress <- aggregate(cbind(breath.count, weight, stress.response,
cort.baseline) ~ week + treatment, stress, mean)
agg_stress_se <- aggregate(cbind(breath.count, weight, stress.response,
cort.baseline) ~ week + treatment, stress, standard_error)
names(agg_stress_se) <- paste(names(agg_stress_se), ".se", sep = "")
agg_stress <- cbind(agg_stress, agg_stress_se[, 3:6])
agg_stress$Week <- 1:4
bs <- 10
gg_breath.count <- ggplot(data = agg_stress, aes(x = Week, y = breath.count,
fill = treatment)) + geom_bar(stat = "identity") + geom_errorbar(aes(ymin = breath.count -
breath.count.se, ymax = breath.count + breath.count.se), width = 0.1) +
scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~treatment,
ncol = 3, scale = "fixed") + labs(y = "Mean change in breath\nrate (breaths/min)",
x = "Week") + theme_classic(base_size = bs) + theme(legend.position = "none")
gg_weight <- ggplot(data = agg_stress, aes(x = Week, y = weight, fill = treatment)) +
geom_bar(stat = "identity") + geom_errorbar(aes(ymin = weight -
weight.se, ymax = weight + weight.se), width = 0.1) + scale_fill_viridis_d(begin = 0.1,
end = 0.9) + facet_wrap(~treatment, ncol = 3, scale = "fixed") +
labs(y = "Mean change in \nweight (grams)", x = "Week") + theme_classic(base_size = bs) +
theme(legend.position = "none")
gg_cort.baseline <- ggplot(data = agg_stress, aes(x = Week, y = cort.baseline,
fill = treatment)) + geom_bar(stat = "identity") + geom_errorbar(aes(ymin = cort.baseline -
cort.baseline.se, ymax = cort.baseline + cort.baseline.se), width = 0.1) +
scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~treatment,
ncol = 3, scale = "fixed") + labs(y = "Mean change in\nbaseline CORT (ng/mL)",
x = "Week") + theme_classic(base_size = bs) + theme(legend.position = "none")
gg_stress.response <- ggplot(data = agg_stress, aes(x = Week, y = stress.response,
fill = treatment)) + geom_bar(stat = "identity") + geom_errorbar(aes(ymin = stress.response -
stress.response.se, ymax = stress.response + stress.response.se),
width = 0.1) + scale_fill_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~treatment, ncol = 3, scale = "fixed") + labs(y = "Mean change in stress response \nCORT (ng/mL)",
x = "Week") + theme_classic(base_size = bs) + theme(legend.position = "none")
gg_coeffs_physio <- lapply(physio_models, function(x) {
vars <- grep("b_", posterior::variables(x), value = TRUE)
draws <- posterior::as_draws_array(x, variable = vars)
coef_table <- draw_summary(draws, variables = vars, probs = c(0.025,
0.975), robust = TRUE)
coef_table$predictor <- rownames(coef_table)
coef_table$predictor <- gsub("b_treatment|b_", "", coef_table$predictor)
coef_table$predictor <- gsub("Stress", " stress", coef_table$predictor)
coef_table$predictor <- gsub("week", "Week", coef_table$predictor)
coef_table <- coef_table[coef_table$predictor != "Intercept",
]
gg_coef <- ggplot2::ggplot(data = coef_table, aes(x = Estimate,
y = predictor)) + geom_vline(xintercept = 0, lty = 2) + ggplot2::geom_point(size = 4,
col = col_pointrange) + ggplot2::geom_errorbar(ggplot2::aes(xmin = `l-95% CI`,
xmax = `u-95% CI`), width = 0, col = col_pointrange) + ggplot2::theme_classic(base_size = bs) +
ggplot2::theme(axis.ticks.length = ggplot2::unit(0, "pt"),
plot.margin = ggplot2::margin(0, 0, 0, 0, "pt"), legend.position = "none",
strip.background = ggplot2::element_blank(), strip.text = ggplot2::element_blank()) +
ggplot2::labs(x = "Effect size", y = "")
return(gg_coef)
})
cowplot::plot_grid(gg_breath.count, gg_weight, gg_cort.baseline, gg_stress.response,
gg_coeffs_physio[[grep("breath", names(gg_coeffs_physio))]] +
theme_classic(base_size = bs), gg_coeffs_physio[[grep("weight",
names(gg_coeffs_physio))]] + theme_classic(base_size = bs),
gg_coeffs_physio[[grep("baseline", names(gg_coeffs_physio))]] +
theme_classic(base_size = bs), gg_coeffs_physio[[grep("response",
names(gg_coeffs_physio))]] + theme_classic(base_size = bs),
nrow = 2, rel_heights = c(1.8, 1))

# try bs = 20 for saving plots
# cowplot::ggsave2(filename =
# './output/bar_graphs_and_estimates_physiology_70dpi.jpeg',
# width = 25, height = 9) cowplot::ggsave2(filename =
# './output/bar_graphs_and_estimates_physiology_300dpi.jpeg',
# dpi = 300, width = 25, height = 9)
Stats
Models: Predicted physio measure ~ treatment + week (continuous) +
IndRandom
Variables (Difference from Week 1): weight, BR, baseline CORT, Stress
CORT, Stress Response
responses <- c("baseline.CORT", "stress.response", "stress.CORT",
"weight", "breath.rate")
predictors <- c("~ treatment + week + (1|ID) + (1|round)")
formulas <- expand.grid(responses = responses, predictors = predictors,
stringsAsFactors = FALSE)
vars_to_scale <- c(responses, "week")
# remove week 1
sub_agg_dat <- agg_dat[agg_dat$week != 1, ]
for (i in vars_to_scale) sub_agg_dat[, vars_to_scale] <- scale(sub_agg_dat[,
vars_to_scale])
physio_models <- lapply(1:nrow(formulas), function(x) {
sub_dat <- sub_agg_dat[!is.na(sub_agg_dat[names(sub_agg_dat) ==
formulas$responses[x]]), ]
sub_dat
mod <- brm(formula = paste(formulas$responses[x], formulas$predictors[x]),
iter = 20000, silent = 2, data = sub_dat, control = list(adapt_delta = 0.9),
chains = 4, prior = c(prior(normal(0, 5), "b"), prior(normal(0,
10), "Intercept"), prior(student_t(3, 0, 10), "sd"), prior(student_t(3,
0, 10), "sigma")))
return(mod)
})
names(physio_models) <- paste(formulas$responses, formulas$predictors)
saveRDS(physio_models, "./data/processed/physiological_response_models.RDS")
physio_models <- readRDS("./data/processed/physiological_response_models.RDS")
for (x in 1:length(physio_models)) html_summary(physio_models[[x]],
gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = FALSE,
remove.intercepts = TRUE, model.name = names(physio_models)[x])
baseline.CORT ~ treatment + week + (1|ID) + (1|round)
|
|
priors
|
formula
|
iterations
|
chains
|
thinning
|
warmup
|
diverg_transitions
|
rhats > 1.05
|
min_bulk_ESS
|
min_tail_ESS
|
seed
|
|
1
|
b-normal(0, 5) Intercept-normal(0, 10) sd-student_t(3, 0, 10)
sigma-student_t(3, 0, 10)
|
baseline.CORT ~ treatment + week + (1 | ID) + (1 | round)
|
20000
|
4
|
1
|
10000
|
406
|
0
|
12756.08
|
17437.15
|
614553873
|
|
|
Estimate
|
l-95% CI
|
u-95% CI
|
Rhat
|
Bulk_ESS
|
Tail_ESS
|
|
HighStress
|
0.760
|
0.108
|
1.420
|
1
|
13169.12
|
17881.30
|
|
MediumStress
|
0.180
|
-0.506
|
0.868
|
1
|
12756.08
|
17437.15
|
|
week
|
-0.001
|
-0.142
|
0.141
|
1
|
27598.67
|
26244.08
|

stress.response ~ treatment + week + (1|ID) + (1|round)
|
|
priors
|
formula
|
iterations
|
chains
|
thinning
|
warmup
|
diverg_transitions
|
rhats > 1.05
|
min_bulk_ESS
|
min_tail_ESS
|
seed
|
|
1
|
b-normal(0, 5) Intercept-normal(0, 10) sd-student_t(3, 0, 10)
sigma-student_t(3, 0, 10)
|
stress.response ~ treatment + week + (1 | ID) + (1 | round)
|
20000
|
4
|
1
|
10000
|
2129
|
0
|
1831.829
|
3145.311
|
1852853265
|
|
|
Estimate
|
l-95% CI
|
u-95% CI
|
Rhat
|
Bulk_ESS
|
Tail_ESS
|
|
HighStress
|
-0.204
|
-0.889
|
0.465
|
1.003
|
1831.829
|
3145.311
|
|
MediumStress
|
0.097
|
-0.667
|
0.839
|
1.002
|
2582.692
|
10967.477
|
|
week
|
-0.152
|
-0.284
|
-0.020
|
1.003
|
6379.514
|
15589.581
|

stress.CORT ~ treatment + week + (1|ID) + (1|round)
|
|
priors
|
formula
|
iterations
|
chains
|
thinning
|
warmup
|
diverg_transitions
|
rhats > 1.05
|
min_bulk_ESS
|
min_tail_ESS
|
seed
|
|
1
|
b-normal(0, 5) Intercept-normal(0, 10) sd-student_t(3, 0, 10)
sigma-student_t(3, 0, 10)
|
stress.CORT ~ treatment + week + (1 | ID) + (1 | round)
|
20000
|
4
|
1
|
10000
|
1015
|
0
|
7696.824
|
8195.241
|
672971911
|
|
|
Estimate
|
l-95% CI
|
u-95% CI
|
Rhat
|
Bulk_ESS
|
Tail_ESS
|
|
HighStress
|
-0.219
|
-0.905
|
0.452
|
1.001
|
8677.295
|
16469.149
|
|
MediumStress
|
0.084
|
-0.684
|
0.832
|
1
|
7696.824
|
8195.241
|
|
week
|
-0.154
|
-0.284
|
-0.022
|
1
|
25064.631
|
25062.987
|

weight ~ treatment + week + (1|ID) + (1|round)
|
|
priors
|
formula
|
iterations
|
chains
|
thinning
|
warmup
|
diverg_transitions
|
rhats > 1.05
|
min_bulk_ESS
|
min_tail_ESS
|
seed
|
|
1
|
b-normal(0, 5) Intercept-normal(0, 10) sd-student_t(3, 0, 10)
sigma-student_t(3, 0, 10)
|
weight ~ treatment + week + (1 | ID) + (1 | round)
|
20000
|
4
|
1
|
10000
|
1142
|
0
|
6459.84
|
11029.72
|
351999546
|
|
|
Estimate
|
l-95% CI
|
u-95% CI
|
Rhat
|
Bulk_ESS
|
Tail_ESS
|
|
HighStress
|
-0.319
|
-0.971
|
0.328
|
1.001
|
10517.09
|
16418.60
|
|
MediumStress
|
-0.120
|
-0.850
|
0.598
|
1.002
|
6459.84
|
11029.72
|
|
week
|
-0.345
|
-0.478
|
-0.211
|
1
|
14142.17
|
16891.28
|

breath.rate ~ treatment + week + (1|ID) + (1|round)
|
|
priors
|
formula
|
iterations
|
chains
|
thinning
|
warmup
|
diverg_transitions
|
rhats > 1.05
|
min_bulk_ESS
|
min_tail_ESS
|
seed
|
|
1
|
b-normal(0, 5) Intercept-normal(0, 10) sd-student_t(3, 0, 10)
sigma-student_t(3, 0, 10)
|
breath.rate ~ treatment + week + (1 | ID) + (1 | round)
|
20000
|
4
|
1
|
10000
|
732
|
0
|
17760.6
|
14209.36
|
1444112821
|
|
|
Estimate
|
l-95% CI
|
u-95% CI
|
Rhat
|
Bulk_ESS
|
Tail_ESS
|
|
HighStress
|
0.147
|
-0.156
|
0.449
|
1
|
17760.60
|
14209.36
|
|
MediumStress
|
0.060
|
-0.264
|
0.389
|
1
|
18526.45
|
24064.77
|
|
week
|
-0.784
|
-0.898
|
-0.670
|
1
|
24048.27
|
26983.78
|

Acoustic space projection
t-SNE
scale_param <- scale(sels[, c("duration", "meanfreq", "sd", "freq.median",
"freq.IQR", "time.IQR", "skew", "kurt", "sp.ent", "time.ent",
"entropy", "meandom", "mindom", "maxdom", "dfrange", "modindx",
"meanpeakf")])
tsne <- Rtsne(scale_param, dims = 2, perplexity = 30, verbose = FALSE,
max_iter = 5000)
saveRDS(tsne, "./data/processed/tsne_on_acoustic_parameters_jun_2021.RDS")
tsne <- readRDS("./data/processed/tsne_on_acoustic_parameters_jun_2021.RDS")
Y <- as.data.frame(tsne$Y)
names(Y) <- c("TSNE1", "TSNE2")
sels <- data.frame(sels, Y)
sels$treatment <- factor(sels$treatment, levels = c("Control", "Medium Stress",
"High Stress"))
ggplot(sels, aes(x = TSNE1, y = TSNE2, col = as.factor(treatment))) +
geom_point() + labs(color = "Treatment") + scale_color_viridis_d(alpha = 0.4) +
theme_classic(base_size = 25) + guides(colour = guide_legend(override.aes = list(size = 10)))

Behavioral parameters
Barplot and effect sizes pgrah
behav_models <- readRDS("./data/processed/behavioral_response_models.RDS")
agg_call.count <- aggregate(cbind(call.count, acoustic.convergence) ~
week + treatment, agg_dat, mean)
agg_behav <- aggregate(cbind(acoustic.diversity, acustic.plasticity) ~
week + treatment, agg_dat, mean)
agg_call.count_se <- aggregate(cbind(call.count, acoustic.convergence) ~
week + treatment, agg_dat, standard_error)
agg_behav_se <- aggregate(cbind(acoustic.diversity, acustic.plasticity) ~
week + treatment, agg_dat, standard_error)
agg_behav_se <- merge(agg_call.count_se, agg_behav_se, all = TRUE)
names(agg_behav_se) <- paste(names(agg_behav_se), ".se", sep = "")
agg_behav <- merge(agg_call.count, agg_behav, all = TRUE)
agg_behav <- cbind(agg_behav, agg_behav_se[, 3:6])
bs <- 10
agg_behav$treatment <- factor(agg_behav$treatment, levels = c("Control",
"Medium Stress", "High Stress"))
gg_call.count <- ggplot(data = agg_behav, aes(x = week, y = call.count,
fill = treatment)) + geom_bar(stat = "identity") + geom_errorbar(aes(ymin = call.count -
call.count.se, ymax = call.count + call.count.se), width = 0.1) +
scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~treatment,
ncol = 3, scale = "fixed") + labs(y = "Vocal output", x = "Week") +
theme_classic(base_size = bs) + theme(legend.position = "none")
gg_acoustic.diversity <- ggplot(data = agg_behav, aes(x = week, y = acoustic.diversity,
fill = treatment)) + geom_bar(stat = "identity") + geom_errorbar(aes(ymin = acoustic.diversity -
acoustic.diversity.se, ymax = acoustic.diversity + acoustic.diversity.se),
width = 0.1) + scale_fill_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~treatment, ncol = 3, scale = "fixed") + labs(y = "Change in vocal diversity",
x = "Week") + theme_classic(base_size = bs) + theme(legend.position = "none")
gg_acustic.plasticity <- ggplot(data = agg_behav, aes(x = week, y = acustic.plasticity,
fill = treatment)) + geom_bar(stat = "identity") + geom_errorbar(aes(ymin = acustic.plasticity -
acustic.plasticity.se, ymax = acustic.plasticity + acustic.plasticity.se),
width = 0.1) + scale_fill_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~treatment, ncol = 3, scale = "fixed") + labs(y = "Vocal plasticity",
x = "Week") + theme_classic(base_size = bs) + theme(legend.position = "none")
gg_acoustic.convergence <- ggplot(data = agg_behav, aes(x = week,
y = acoustic.convergence, fill = treatment)) + geom_bar(stat = "identity") +
geom_errorbar(aes(ymin = acoustic.convergence - acoustic.convergence.se,
ymax = acoustic.convergence + acoustic.convergence.se), width = 0.1) +
scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~treatment,
ncol = 3, scale = "fixed") + labs(y = "Vocal convergence", x = "Week") +
theme_classic(base_size = bs) + theme(legend.position = "none")
gg_coeffs_behav <- lapply(behav_models, function(x) {
vars <- grep("b_", posterior::variables(x), value = TRUE)
draws <- posterior::as_draws_array(x, variable = vars)
coef_table <- draw_summary(draws, variables = vars, probs = c(0.025,
0.975), robust = TRUE)
coef_table$predictor <- rownames(coef_table)
coef_table$predictor <- gsub("b_treatment|b_", "", coef_table$predictor)
coef_table$predictor <- gsub("Stress", " stress", coef_table$predictor)
coef_table$predictor <- gsub("week", "Week", coef_table$predictor)
coef_table <- coef_table[coef_table$predictor != "Intercept",
]
gg_coef <- ggplot2::ggplot(data = coef_table, aes(x = Estimate,
y = predictor)) + geom_vline(xintercept = 0, lty = 2) + ggplot2::geom_point(size = 4,
col = col_pointrange) + ggplot2::geom_errorbar(ggplot2::aes(xmin = `l-95% CI`,
xmax = `u-95% CI`), width = 0, col = col_pointrange) + ggplot2::theme_classic(base_size = bs) +
ggplot2::theme(axis.ticks.length = ggplot2::unit(0, "pt"),
plot.margin = ggplot2::margin(0, 0, 0, 0, "pt"), legend.position = "none",
strip.background = ggplot2::element_blank(), strip.text = ggplot2::element_blank()) +
ggplot2::labs(x = "Effect size", y = "")
return(gg_coef)
})
cowplot::plot_grid(gg_call.count, gg_acoustic.diversity, gg_acustic.plasticity,
gg_acoustic.convergence, gg_coeffs_behav[[grep("count", names(gg_coeffs_behav))]] +
theme_classic(base_size = bs), gg_coeffs_behav[[grep("diversity",
names(gg_coeffs_behav))]] + theme_classic(base_size = bs),
gg_coeffs_behav[[grep("plasticity", names(gg_coeffs_behav))]] +
theme_classic(base_size = bs), gg_coeffs_behav[[grep("convergence",
names(gg_coeffs_behav))]] + theme_classic(base_size = bs),
nrow = 2, rel_heights = c(1.8, 1))

# try bs = 20 for saving plots
# cowplot::ggsave2(filename =
# './output/bar_graphs_and_estimates_behavior_70dpi.jpeg', width
# = 25, height = 9) # cowplot::ggsave2(filename =
# './output/bar_graphs_and_estimates_behavior_300dpi.jpeg', dpi
# = 300, width = 25, height = 9)
Stats
Model: Predicted behavior ~ treatment + week (continuous) +
IndRandom
Variables: # calls, Distance moved from self in first week, Overlap
to original acoustic space, Match to group repertoire, Maybe overall
size of acoustic space
Do as comparison to week one using rarefacted calls and kernel
density
responses <- c("call.count", "acoustic.diversity", "acustic.plasticity",
"acoustic.convergence")
predictors <- c("~ treatment + week + (1|ID) + (1|round)")
formulas <- expand.grid(responses = responses, predictors = predictors,
stringsAsFactors = FALSE)
vars_to_scale <- c(responses, "week")
for (i in vars_to_scale) agg_dat[, vars_to_scale] <- scale(agg_dat[,
vars_to_scale])
behav_models <- lapply(1:nrow(formulas), function(x) {
sub_dat <- agg_dat[!is.na(agg_dat[names(agg_dat) == formulas$responses[x]]),
]
# remove week 1
if (!grepl("count|group", formulas$responses[x]))
sub_dat <- sub_dat[sub_dat$week != 1, ]
mod <- brm(formula = paste(formulas$responses[x], formulas$predictors[x]),
iter = 50000, silent = 2, data = sub_dat, control = list(adapt_delta = 0.9,
max_treedepth = 15), chains = 4, prior = c(prior(normal(0,
5), "b"), prior(normal(0, 10), "Intercept"), prior(student_t(3,
0, 10), "sd"), prior(student_t(3, 0, 10), "sigma")))
return(mod)
})
names(behav_models) <- paste(formulas$responses, formulas$predictors)
saveRDS(behav_models, "./data/processed/behavioral_response_models.RDS")
behav_models <- readRDS("./data/processed/behavioral_response_models.RDS")
for (x in 1:length(behav_models)) html_summary(behav_models[[x]],
gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = FALSE,
remove.intercepts = TRUE, model.name = names(behav_models)[x])
call.count ~ treatment + week + (1|ID) + (1|round)
|
|
priors
|
formula
|
iterations
|
chains
|
thinning
|
warmup
|
diverg_transitions
|
rhats > 1.05
|
min_bulk_ESS
|
min_tail_ESS
|
seed
|
|
1
|
b-normal(0, 5) Intercept-normal(0, 10) sd-student_t(3, 0, 10)
sigma-student_t(3, 0, 10)
|
call.count ~ treatment + week + (1 | ID) + (1 | round)
|
50000
|
4
|
1
|
25000
|
1455
|
0
|
21916.29
|
19215.21
|
2086259884
|
|
|
Estimate
|
l-95% CI
|
u-95% CI
|
Rhat
|
Bulk_ESS
|
Tail_ESS
|
|
HighStress
|
0.784
|
0.205
|
1.357
|
1
|
21916.29
|
19215.21
|
|
MediumStress
|
0.394
|
-0.183
|
0.971
|
1
|
29122.91
|
29122.87
|
|
week
|
-0.059
|
-0.203
|
0.086
|
1
|
48490.61
|
55075.00
|

acoustic.diversity ~ treatment + week + (1|ID) + (1|round)
|
|
priors
|
formula
|
iterations
|
chains
|
thinning
|
warmup
|
diverg_transitions
|
rhats > 1.05
|
min_bulk_ESS
|
min_tail_ESS
|
seed
|
|
1
|
b-normal(0, 5) Intercept-normal(0, 10) sd-student_t(3, 0, 10)
sigma-student_t(3, 0, 10)
|
acoustic.diversity ~ treatment + week + (1 | ID) + (1 | round)
|
50000
|
4
|
1
|
25000
|
2953
|
0
|
6138.296
|
3258.34
|
1000749341
|
|
|
Estimate
|
l-95% CI
|
u-95% CI
|
Rhat
|
Bulk_ESS
|
Tail_ESS
|
|
HighStress
|
-0.208
|
-1.070
|
0.657
|
1.001
|
6138.296
|
3258.34
|
|
MediumStress
|
-0.699
|
-1.559
|
0.158
|
1.001
|
8701.168
|
12053.83
|
|
week
|
0.036
|
-0.081
|
0.152
|
1
|
27120.995
|
32683.11
|

acustic.plasticity ~ treatment + week + (1|ID) + (1|round)
|
|
priors
|
formula
|
iterations
|
chains
|
thinning
|
warmup
|
diverg_transitions
|
rhats > 1.05
|
min_bulk_ESS
|
min_tail_ESS
|
seed
|
|
1
|
b-normal(0, 5) Intercept-normal(0, 10) sd-student_t(3, 0, 10)
sigma-student_t(3, 0, 10)
|
acustic.plasticity ~ treatment + week + (1 | ID) + (1 | round)
|
50000
|
4
|
1
|
25000
|
7397
|
0
|
1025.75
|
483.942
|
803966074
|
|
|
Estimate
|
l-95% CI
|
u-95% CI
|
Rhat
|
Bulk_ESS
|
Tail_ESS
|
|
HighStress
|
-0.904
|
-1.731
|
-0.062
|
1.003
|
4966.382
|
30252.390
|
|
MediumStress
|
-0.232
|
-1.130
|
0.678
|
1.005
|
1025.750
|
483.942
|
|
week
|
0.171
|
0.034
|
0.306
|
1.002
|
2427.679
|
26782.955
|

acoustic.convergence ~ treatment + week + (1|ID) + (1|round)
|
|
priors
|
formula
|
iterations
|
chains
|
thinning
|
warmup
|
diverg_transitions
|
rhats > 1.05
|
min_bulk_ESS
|
min_tail_ESS
|
seed
|
|
1
|
b-normal(0, 5) Intercept-normal(0, 10) sd-student_t(3, 0, 10)
sigma-student_t(3, 0, 10)
|
acoustic.convergence ~ treatment + week + (1 | ID) + (1 | round)
|
50000
|
4
|
1
|
25000
|
8224
|
0
|
150.4
|
57.174
|
1972081133
|
|
|
Estimate
|
l-95% CI
|
u-95% CI
|
Rhat
|
Bulk_ESS
|
Tail_ESS
|
|
HighStress
|
0.335
|
-0.307
|
1.003
|
1.006
|
1086.629
|
21463.246
|
|
MediumStress
|
0.372
|
-0.288
|
1.045
|
1.003
|
12829.194
|
16680.338
|
|
week
|
-0.228
|
-0.406
|
-0.058
|
1.017
|
150.400
|
57.174
|

Takeaways
Higher vocal output in “high stress” birds compared to
control
Higher acoustic overlap to themselves in week 1 for “high stress”
birds compared to control
Decrease in overlap to themselves in week 1 across time
Combined model diagnostics
check_rds_models(path = "./data/processed", html = TRUE, verbose = FALSE)
|
|
priors
|
formula
|
iterations
|
chains
|
thinning
|
warmup
|
diverg_transitions
|
rhats > 1.05
|
min_bulk_ESS
|
min_tail_ESS
|
seed
|
|
12
|
Intercept-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5)
|
size ~ treatment + stress.response
|
5000
|
1
|
1
|
2500
|
0
|
0
|
1055.507
|
1103.504
|
654213680
|
|
13
|
Intercept-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5)
|
size ~ treatment + stress.response
|
5000
|
1
|
1
|
2500
|
0
|
0
|
807.903
|
988.693
|
5788341
|
|
14
|
Intercept-student_t(3, 422, 394.4) sigma-student_t(3, 0, 394.4)
|
selec ~ treatment + stress.response
|
5000
|
1
|
1
|
2500
|
0
|
0
|
439.266
|
695.445
|
825114723
|
|
15
|
Intercept-student_t(3, 111, 143.8) sigma-student_t(3, 0, 143.8)
|
selec ~ treatment + stress.response
|
5000
|
1
|
1
|
2500
|
0
|
0
|
641.073
|
841.561
|
1461380493
|
|
16
|
Intercept-student_t(3, 0, 2.5) phi-gamma(0.01, 0.01)
|
overlap.to.first.week ~ treatment + stress.response
|
5000
|
1
|
1
|
2500
|
0
|
0
|
1133.767
|
1203.772
|
363627053
|
|
17
|
Intercept-student_t(3, 0, 2.5) phi-gamma(0.01, 0.01)
|
overlap.to.first.week ~ treatment + stress.response
|
5000
|
1
|
1
|
2500
|
0
|
0
|
1196.878
|
1623.758
|
1435724951
|
Session information
## R version 4.1.0 (2021-05-18)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
## LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=pt_BR.UTF-8
## [7] LC_PAPER=es_CR.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] PhenotypeSpace_0.1.0 warbleR_1.1.26 NatureSounds_1.0.4
## [4] seewave_2.1.8 tuneR_1.3.3.1 ggridges_0.5.3
## [7] posterior_1.2.0 ggrepel_0.9.1 cowplot_1.1.1
## [10] tidybayes_3.0.2 modelr_0.1.8 tidyr_1.2.0
## [13] forcats_0.5.1 purrr_0.3.4 dplyr_1.0.8
## [16] lme4_1.1-28 Matrix_1.3-4 brms_2.16.3
## [19] Rcpp_1.0.8 GGally_2.1.2 sp_1.4-6
## [22] MASS_7.3-54 formatR_1.11 knitr_1.37
## [25] kableExtra_1.3.4 ggplot2_3.3.5 viridis_0.6.2
## [28] viridisLite_0.4.0 pbapply_1.5-0 readxl_1.3.1
## [31] lubridate_1.8.0 remotes_2.4.2
##
## loaded via a namespace (and not attached):
## [1] utf8_1.2.2 tidyselect_1.1.2 fftw_1.0-6.1
## [4] htmlwidgets_1.5.4 grid_4.1.0 munsell_0.5.0
## [7] codetools_0.2-18 DT_0.20 miniUI_0.1.1.1
## [10] withr_2.4.3 spatstat.random_2.1-0 Brobdingnag_1.2-7
## [13] colorspace_2.0-3 highr_0.9 uuid_1.1-0
## [16] rstudioapi_0.13 stats4_4.1.0 dtw_1.22-3
## [19] tensor_1.5 bayesplot_1.8.1 labeling_0.4.2
## [22] emmeans_1.8.1-1 rstan_2.21.3 polyclip_1.10-0
## [25] farver_2.1.0 bridgesampling_1.1-2 rprojroot_2.0.2
## [28] coda_0.19-4 vctrs_0.3.8 generics_0.1.2
## [31] TH.data_1.1-0 xfun_0.29 R6_2.5.1
## [34] markdown_1.1 gamm4_0.2-6 projpred_2.0.2
## [37] bitops_1.0-7 spatstat.utils_2.3-0 reshape_0.8.8
## [40] assertthat_0.2.1 promises_1.2.0.1 scales_1.1.1
## [43] multcomp_1.4-17 gtable_0.3.0 goftest_1.2-3
## [46] processx_3.5.2 xaringanExtra_0.7.0 sandwich_3.0-1
## [49] rlang_1.0.1 systemfonts_1.0.4 splines_4.1.0
## [52] spatstat.geom_2.3-2 broom_0.7.12 checkmate_2.0.0
## [55] inline_0.3.19 yaml_2.3.5 reshape2_1.4.4
## [58] abind_1.4-5 threejs_0.3.3 crosstalk_1.2.0
## [61] backports_1.4.1 httpuv_1.6.5 rsconnect_0.8.25
## [64] tensorA_0.36.2 tools_4.1.0 ellipsis_0.3.2
## [67] spatstat.core_2.4-0 raster_3.5-15 jquerylib_0.1.4
## [70] RColorBrewer_1.1-2 proxy_0.4-26 plyr_1.8.6
## [73] base64enc_0.1-3 RCurl_1.98-1.6 ps_1.6.0
## [76] prettyunits_1.1.1 rpart_4.1-15 deldir_1.0-6
## [79] zoo_1.8-9 cluster_2.1.2 magrittr_2.0.2
## [82] ggdist_3.1.0 colourpicker_1.1.1 mvtnorm_1.1-3
## [85] matrixStats_0.61.0 shinyjs_2.1.0 mime_0.12
## [88] evaluate_0.15 arrayhelpers_1.1-0 xtable_1.8-4
## [91] shinystan_2.5.0 gridExtra_2.3 rstantools_2.1.1
## [94] compiler_4.1.0 tibble_3.1.6 crayon_1.5.0
## [97] minqa_1.2.4 StanHeaders_2.21.0-7 htmltools_0.5.2
## [100] mgcv_1.8-36 later_1.3.0 RcppParallel_5.1.5
## [103] DBI_1.1.1 boot_1.3-28 permute_0.9-7
## [106] cli_3.2.0 parallel_4.1.0 igraph_1.2.11
## [109] pkgconfig_2.0.3 signal_0.7-7 terra_1.5-21
## [112] spatstat.sparse_2.1-0 xml2_1.3.3 svUnit_1.0.6
## [115] dygraphs_1.1.1.6 svglite_2.1.0 bslib_0.3.1
## [118] webshot_0.5.2 estimability_1.4.1 rvest_1.0.2
## [121] stringr_1.4.0 distributional_0.3.0 callr_3.7.0
## [124] digest_0.6.29 vegan_2.5-7 spatstat.data_2.1-2
## [127] rmarkdown_2.11 cellranger_1.1.0 shiny_1.7.1
## [130] gtools_3.9.2 rjson_0.2.21 nloptr_2.0.0
## [133] lifecycle_1.0.1 nlme_3.1-152 jsonlite_1.8.0
## [136] fansi_1.0.2 pillar_1.7.0 lattice_0.20-44
## [139] loo_2.4.1 fastmap_1.1.0 httr_1.4.2
## [142] pkgbuild_1.3.1 survival_3.2-11 glue_1.6.1
## [145] xts_0.12.1 shinythemes_1.2.0 stringi_1.7.6
## [148] sass_0.4.0